% This script plots the simulated social welfare function described
% in "Optimal Defaults with Normative Ambiguity" by Jacob Goldin and Daniel Reck.

clear
clc

%cd into the directory where you want to run the code
%cd 'XXX/finalcode'

%Parameters

%Preferences
alpha=.5;

%Distributions
sigma=1;
xm=0;
gamma=1;

%incidental
xi=sqrt(2*gamma/alpha);

%Range of default
d=-5:0.1:5;

%setup loop
W=zeros(1,length(d));
dW=zeros(1,length(d));



%pi=1
pi=1;
%evaluate welfare
for i=1:length(d)
    thisd=d(i);
    xl=thisd-xi;
    xu=thisd+xi;
    W(i)=   integral(@(x) -pi*gamma*normpdf(x,xm,sigma),-100000, xl) + ...
            integral(@(x) -alpha/2*(thisd-x).^2.*normpdf(x,xm,sigma),xl,xu)+ ...
            integral(@(x) -pi*gamma*normpdf(x,xm,sigma),xu,100000);
    dW(i)=integral(@(x) -alpha/2*(thisd-x).*normpdf(x,xm,sigma),xl,xu) + (1-pi)*gamma*(normpdf(xl,xm,sigma)-normpdf(xu,xm,sigma));
end

plot(d,W)

figure('Color',[1 1 1]);
hold on
set(gca, 'fontsize',14);
plot(d, W, 'LineWidth',2)
xlabel('d');
xticks(-5:1:5)
ylabel('W(d)');
savefig('output/Wd_pi1')
print('output/Wd_pi1','-djpeg')
hold off


%pi=0
pi=0;
%evaluate welfare
for i=1:length(d)
    thisd=d(i);
    xl=thisd-xi;
    xu=thisd+xi;
    W(i)=   integral(@(x) -pi*gamma*normpdf(x,xm,sigma),-100000, xl) + ...
            integral(@(x) -alpha/2*(thisd-x).^2.*normpdf(x,xm,sigma),xl,xu)+ ...
            integral(@(x) -pi*gamma*normpdf(x,xm,sigma),xu,100000);
    dW(i)=integral(@(x) -alpha/2*(thisd-x).*normpdf(x,xm,sigma),xl,xu) + (1-pi)*gamma*(normpdf(xl,xm,sigma)-normpdf(xu,xm,sigma));
end

plot(d,W)

figure('Color',[1 1 1]);
hold on
set(gca, 'fontsize',14);
plot(d, W, 'LineWidth',2)
xlabel('d');
xticks(-5:1:5)
ylabel('W(d)');
savefig('output/Wd_pi0')
print('output/Wd_pi0','-djpeg')
hold off

%pi=pibar=.19;
pi=0.19;
%evaluate welfare
for i=1:length(d)
    thisd=d(i);
    xl=thisd-xi;
    xu=thisd+xi;
    W(i)=   integral(@(x) -pi*gamma*normpdf(x,xm,sigma),-100000, xl) + ...
            integral(@(x) -alpha/2*(thisd-x).^2.*normpdf(x,xm,sigma),xl,xu)+ ...
            integral(@(x) -pi*gamma*normpdf(x,xm,sigma),xu,100000);
    dW(i)=integral(@(x) -alpha/2*(thisd-x).*normpdf(x,xm,sigma),xl,xu) + (1-pi)*gamma*(normpdf(xl,xm,sigma)-normpdf(xu,xm,sigma));
end

plot(d,W)

figure('Color',[1 1 1]);
hold on
set(gca, 'fontsize',14);
plot(d, W, 'LineWidth',2)
xlabel('d');
xticks(-5:1:5)
ylabel('W(d)');
savefig('output/Wd_pibar')
print('output/Wd_pibar','-djpeg')
hold off

%Loop over pi

piv=0:0.2:1;
W=zeros(length(d),length(piv));
for j=1:length(piv)
for i=1:length(d)
    pi=piv(j);
    thisd=d(i);
    xl=thisd-xi;
    xu=thisd+xi;
    W(i,j)=   integral(@(x) -pi*gamma*normpdf(x,xm,sigma),-100000, xl) + ...
            integral(@(x) -alpha/2*(thisd-x).^2.*normpdf(x,xm,sigma),xl,xu)+ ...
            integral(@(x) -pi*gamma*normpdf(x,xm,sigma),xu,100000);
end
end

plot(d, W(1:length(d),:),'LineWidth',2)
legend('pi=0','pi=0.2','pi=0.4','pi=0.6','pi=0.8','pi=1.0',...
    'Location','southoutside','Orientation','horozontal');
xlabel('d');
ylabel('W(d)');
savefig('output/fig1_Wd_bypi')
print('output/fig1_Wd_bypi','-djpeg')
hold off

